---
title: "Placebo thresholds"
output: pdf_document
date: "2025-01-13"
---

```{r , setup, include=FALSE, eval = FALSE}
knitr::opts_chunk$set(
  comment = '', warning=FALSE,  message=FALSE, results= FALSE , fig.width = 9,fig.height = 7, echo=FALSE, root.dir =""
  
)

knitr::opts_knit$set(
  root.dir =""
)

```

```{r packages}
rm(list=ls())
  library(tidyverse)
  library(gridExtra)
  library(dplyr)
  library(rdrobust)
  library(rddensity)
  library(viridisLite)
  library(foreign)
  library(ggplot2)
  library(data.table)
  library(doBy)
  library(socviz)
  library(openxlsx) 
  library(modelsummary)
  library(labelled)
  library(gtsummary)
  library(gt)
  library(kableExtra)
  library(ggpubr)
```

```{r loaddata}
data_treatment<-data.frame(read.csv("1.data/2.intermediatedata/carrs_full.csv"))
data_treatment2<-data_treatment[data_treatment$female==0,]
data_treatment1<-data_treatment[data_treatment$female==1,]


data_diagnosis<-data.frame(read.csv("1.data/2.intermediatedata/carrs_selected.csv"))
data_diagnosis2<-data_diagnosis[data_diagnosis$female==0,]
data_diagnosis1<-data_diagnosis[data_diagnosis$female==1,]

data_bp<-data.frame(read.csv("1.data/2.intermediatedata/carrs_bp.csv"))
data_bp2<-data_bp[data_bp$female==0,]
data_bp1<-data_bp[data_bp$female==1,]

# lists with placebo thresholds (-0.1, -0.05, 0, 0.05, 0.1)
  thresh = rep(seq(from=-0.1, to=0.1, by=0.05), 1)
  thresh
```

``` {r PTcat, eval = TRUE}
 pthresh <- function(Yi,pt) {
    out =  rdrobust(y = Yi, X, kernel = "triangular", p = 1, all = T,
                    c = pt, cluster = clus, covs = cov, rho = 1 )
    summary(out)
    
    ret <- data.frame(
      res = c(as.character(round(100*out$coef[1,1],2)), as.character(round(out$pv[3,1],2)), 
              paste0( "(", as.character(round(100*out$ci[3,1],2)), " - ", as.character(round(100*out$ci[3,2],2)), ")" ), 
              as.character(round(out$bws[1,1],2)), as.character(out$N_h[1] + out$N_h[2]))
    )
    rownames(ret) <- c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")
    ret
  
    return(list(ret=ret))
  }
```

``` {r PTcont, eval = TRUE}
 pthresh_cont <- function(Yi,pt) {
    out =  rdrobust(y = Yi, X, kernel = "triangular", p = 1, all = T,
                    c = pt, cluster = clus, covs = cov, rho = 1 )
    summary(out)
    
    ret <- data.frame(
      res = c(as.character(round(out$coef[1,1],2)), as.character(round(out$pv[3,1],2)), 
              paste0( "(", as.character(round(out$ci[3,1],2)), " - ", as.character(round(out$ci[3,2],2)), ")" ), 
              as.character(round(out$bws[1,1],2)), as.character(out$N_h[1] + out$N_h[2]))
    )
    rownames(ret) <- c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")
    ret
  
    return(list(ret=ret))
  }
```

``` {r PlaceboThresholdsDiag, eval = TRUE }

## Full sample
  X = data_diagnosis$runningvar
  Y = data_diagnosis$hypdrug
  clus = data_diagnosis$pid
  cov = cbind(data_diagnosis$w01, data_diagnosis$w23, data_diagnosis$w45)
  
  PTDT <- data.frame(matrix(, nrow=5, ncol=0))
  rownames(PTDT) <- c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")
  for (i in 1:length(thresh)) {
    assign("PTDTres", pthresh(Yi = Y,pt=thresh[i]))
    PTDT[i] <- PTDTres
  }
  
## Men
  X = data_diagnosis2$runningvar_male
  Y = data_diagnosis2$hypdrug
  clus = data_diagnosis2$pid
  cov = cbind(data_diagnosis2$w01, data_diagnosis2$w23, data_diagnosis2$w45)
  
  PTDM <- data.frame(matrix(, nrow=5, ncol=0))
  rownames(PTDM) <- c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")
  for (i in 1:length(thresh)) {
    assign("PTDMres", pthresh(Yi = Y,pt=thresh[i]))
    PTDM[i] <- PTDMres
  }

## Women
  X = data_diagnosis1$runningvar_female
  Y = data_diagnosis1$hypdrug
  clus = data_diagnosis1$pid
  cov = cbind(data_diagnosis1$w01, data_diagnosis1$w23, data_diagnosis1$w45)
  
    PTDF <- data.frame(matrix(, nrow=5, ncol=0))
  rownames(PTDF) <- c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")
  for (i in 1:length(thresh)) {
    assign("PTDFres", pthresh(Yi = Y,pt=thresh[i]))
    PTDF[i] <- PTDFres
  }

``` 


``` {r PlaceboThresholdsTreat, eval = TRUE }
## Full sample
  X = data_treatment$runningvar
  Y = data_treatment$hypdrug
  clus = data_treatment$pid
  cov = cbind(data_treatment$w01, data_treatment$w23, data_treatment$w45)
  
  PTTT <- data.frame(matrix(, nrow=5, ncol=0))
  rownames(PTTT) <- c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")
  for (i in 1:length(thresh)) {
    assign("PTTTres", pthresh(Yi = Y,pt=thresh[i]))
    PTTT[i] <- PTTTres
  }
  
## Men
  X = data_treatment2$runningvar_male
  Y = data_treatment2$hypdrug
  clus = data_treatment2$pid
  cov = cbind(data_treatment2$w01, data_treatment2$w23, data_treatment2$w45)
  
  PTTM <- data.frame(matrix(, nrow=5, ncol=0))
  rownames(PTTM) <- c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")
  for (i in 1:length(thresh)) {
    assign("PTTMres", pthresh(Yi = Y,pt=thresh[i]))
    PTTM[i] <- PTTMres
  }

## Women
  X = data_treatment1$runningvar_female
  Y = data_treatment1$hypdrug
  clus = data_treatment1$pid
  cov = cbind(data_treatment1$w01, data_treatment1$w23, data_treatment1$w45)
  
    PTTF <- data.frame(matrix(, nrow=5, ncol=0))
  rownames(PTTF) <- c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")
  for (i in 1:length(thresh)) {
    assign("PTTFres", pthresh(Yi = Y,pt=thresh[i]))
    PTTF[i] <- PTTFres
  }

``` 

``` {r PlaceboThresholdsSBP, eval = TRUE }
## Full sample
  X = data_bp$runningvar
  Y = data_bp$sysdiff
  clus = data_bp$pid
  cov = cbind(data_bp$w02, data_bp$w24)
  
  PTST <- data.frame(matrix(, nrow=5, ncol=0))
  rownames(PTST) <- c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")
  for (i in 1:length(thresh)) {
    assign("PTSTres", pthresh_cont(Yi = Y,pt=thresh[i]))
    PTST[i] <- PTSTres
  }
  
## Men
  X = data_bp2$runningvar_male
  Y = data_bp2$sysdiff
  clus = data_bp2$pid
  cov = cbind(data_bp2$w02, data_bp2$w24)
  
  PTSM <- data.frame(matrix(, nrow=5, ncol=0))
  rownames(PTSM) <- c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")
  for (i in 1:length(thresh)) {
    assign("PTSMres", pthresh_cont(Yi = Y,pt=thresh[i]))
    PTSM[i] <- PTSMres
  }

## Women
  X = data_bp1$runningvar_female
  Y = data_bp1$sysdiff
  clus = data_bp1$pid
  cov = cbind(data_bp1$w02, data_bp1$w24)
  
  PTSF <- data.frame(matrix(, nrow=5, ncol=0))
  rownames(PTSF) <- c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")
  for (i in 1:length(thresh)) {
    assign("PTSFres", pthresh_cont(Yi = Y,pt=thresh[i]))
    PTSF[i] <- PTSFres
  }

``` 

``` {r PlaceboThresholdsDBP, eval = TRUE }
## Full sample
  X = data_bp$runningvar
  Y = data_bp$diasdiff
  clus = data_bp$pid
  cov = cbind(data_bp$w02, data_bp$w24)
  
  PTDT <- data.frame(matrix(, nrow=5, ncol=0))
  rownames(PTDT) <- c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")
  for (i in 1:length(thresh)) {
    assign("PTDTres", pthresh_cont(Yi = Y,pt=thresh[i]))
    PTDT[i] <- PTDTres
  }
  
## Men
  X = data_bp2$runningvar_male
  Y = data_bp2$diasdiff
  clus = data_bp2$pid
  cov = cbind(data_bp2$w02, data_bp2$w24)
  
  PTDM <- data.frame(matrix(, nrow=5, ncol=0))
  rownames(PTDM) <- c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")
  for (i in 1:length(thresh)) {
    assign("PTDMres", pthresh_cont(Yi = Y,pt=thresh[i]))
    PTDM[i] <- PTDMres
  }

## Women
  X = data_bp1$runningvar_female
  Y = data_bp1$diasdiff
  clus = data_bp1$pid
  cov = cbind(data_bp1$w02, data_bp1$w24)
  
  PTDF <- data.frame(matrix(, nrow=5, ncol=0))
  rownames(PTDF) <- c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")
  for (i in 1:length(thresh)) {
    assign("PTDFres", pthresh_cont(Yi = Y,pt=thresh[i]))
    PTDF[i] <- PTDFres
  }

``` 

```{r tablePTD}
##COMBINED - DIAGNOSIS
  filename <- paste0("pt_diag_combined.tex")
  caption=paste0("Effect of home-based screening on hypertension diagnosis with placebo thresholds")

comb_PTD<-rbind(as.matrix(PTDT),as.matrix(PTDF),as.matrix(PTDM))

combined<-kbl(comb_PTD,caption=caption, booktabs = T, linesep = '', escape = FALSE,align = "lccccc", format = "latex",digits = 2, col.names = c("(1)", "(2)", "(3)", "(4)","(5)")) %>%
    kable_styling(latex_options = "HOLD_position", position = "center") %>%
    pack_rows("Panel A: Total", 1,5) %>%
    pack_rows("Estimation results", 1,3) %>%
    pack_rows("Bandwidth details", 4, 5) %>%
    pack_rows(" ", 6,10) %>%
    pack_rows("Panel B: Female", 6,10) %>%
    pack_rows("Estimation results", 6, 8) %>%
    pack_rows("Bandwidth details", 9, 10)  %>%
    pack_rows(" ", 11,15) %>%
    pack_rows("Panel C: Male", 11,15) %>%
    pack_rows("Estimation results", 11, 13) %>%
    pack_rows("Bandwidth details", 14, 15)  %>%
    footnote(general = c("The model estimation included a local linear trend, triangular Kernel weights, and a data-driven bandwidth selection (mean squared error optimal bandwidth). The table displays the local average treatment effect in percentage points for individuals with a standardized blood pressure of zero and CIs resulting from robust bias-corrected standard errors clustered at the individual-level. The bandwidth was varied by +/- 0.2 in intervals of 0.05, column (5) contains the results."," Abbreviation: CI = confidence interval"),general_title="",threeparttable = T, fixed_small_size = T)%>%
    save_kable(filename, format = "latex")

```

```{r tablePTT}
##COMBINED - Treatment
  filename <- paste0("pt_treat_combined.tex")
  caption=paste0("Effect of home-based screening on hypertension treatment with placebo thresholds")

comb_PTT<-rbind(as.matrix(PTTT),as.matrix(PTTF),as.matrix(PTTM))

combined<-kbl(comb_PTT,caption=caption, booktabs = T, linesep = '', escape = FALSE,align = "lccccc", format = "latex",digits = 2, col.names = c("(1)", "(2)", "(3)", "(4)","(5)")) %>%
    kable_styling(latex_options = "HOLD_position", position = "center") %>%
    pack_rows("Panel A: Total", 1,5) %>%
    pack_rows("Estimation results", 1,3) %>%
    pack_rows("Bandwidth details", 4, 5) %>%
    pack_rows(" ", 6,10) %>%
    pack_rows("Panel B: Female", 6,10) %>%
    pack_rows("Estimation results", 6, 8) %>%
    pack_rows("Bandwidth details", 9, 10)  %>%
    pack_rows(" ", 11,15) %>%
    pack_rows("Panel C: Male", 11,15) %>%
    pack_rows("Estimation results", 11, 13) %>%
    pack_rows("Bandwidth details", 14, 15)  %>%
    footnote(general = c("The model estimation included a local linear trend, triangular Kernel weights, and a data-driven bandwidth selection (mean squared error optimal bandwidth). The table displays the local average treatment effect in percentage points for individuals with a standardized blood pressure of zero and CIs resulting from robust bias-corrected standard errors clustered at the individual-level. The bandwidth was varied by +/- 0.2 in intervals of 0.05, column (5) contains the results."," Abbreviation: CI = confidence interval"),general_title="",threeparttable = T, fixed_small_size = T)%>%
    save_kable(filename, format = "latex")

```

```{r tablePTS}
##COMBINED - Treatment
  filename <- paste0("pt_sbp_combined.tex")
  caption=paste0("Effect of home-based screening on systolic blood pressure with placebo thresholds")

comb_PTS<-rbind(as.matrix(PTST),as.matrix(PTSF),as.matrix(PTSM))

combined<-kbl(comb_PTS,caption=caption, booktabs = T, linesep = '', escape = FALSE,align = "lccccc", format = "latex",digits = 2, col.names = c("(1)", "(2)", "(3)", "(4)","(5)")) %>%
    kable_styling(latex_options = "HOLD_position", position = "center") %>%
    pack_rows("Panel A: Total", 1,5) %>%
    pack_rows("Estimation results", 1,3) %>%
    pack_rows("Bandwidth details", 4, 5) %>%
    pack_rows(" ", 6,10) %>%
    pack_rows("Panel B: Female", 6,10) %>%
    pack_rows("Estimation results", 6, 8) %>%
    pack_rows("Bandwidth details", 9, 10)  %>%
    pack_rows(" ", 11,15) %>%
    pack_rows("Panel C: Male", 11,15) %>%
    pack_rows("Estimation results", 11, 13) %>%
    pack_rows("Bandwidth details", 14, 15)  %>%
    footnote(general = c("The model estimation included a local linear trend, triangular Kernel weights, and a data-driven bandwidth selection (mean squared error optimal bandwidth). The table displays the local average treatment effect in percentage points for individuals with a standardized blood pressure of zero and CIs resulting from robust bias-corrected standard errors clustered at the individual-level. The bandwidth was varied by +/- 0.2 in intervals of 0.05, column (5) contains the results."," Abbreviation: CI = confidence interval"),general_title="",threeparttable = T, fixed_small_size = T)%>%
    save_kable(filename, format = "latex")

```

```{r tablePTD}
##COMBINED - Treatment
  filename <- paste0("pt_dbp_combined.tex")
  caption=paste0("Effect of home-based screening on diastolic blood pressure with placebo thresholds")

comb_PTD<-rbind(as.matrix(PTDT),as.matrix(PTDF),as.matrix(PTDM))

combined<-kbl(comb_PTD,caption=caption, booktabs = T, linesep = '', escape = FALSE,align = "lccccc", format = "latex",digits = 2, col.names = c("(1)", "(2)", "(3)", "(4)","(5)")) %>%
    kable_styling(latex_options = "HOLD_position", position = "center") %>%
    pack_rows("Panel A: Total", 1,5) %>%
    pack_rows("Estimation results", 1,3) %>%
    pack_rows("Bandwidth details", 4, 5) %>%
    pack_rows(" ", 6,10) %>%
    pack_rows("Panel B: Female", 6,10) %>%
    pack_rows("Estimation results", 6, 8) %>%
    pack_rows("Bandwidth details", 9, 10)  %>%
    pack_rows(" ", 11,15) %>%
    pack_rows("Panel C: Male", 11,15) %>%
    pack_rows("Estimation results", 11, 13) %>%
    pack_rows("Bandwidth details", 14, 15)  %>%
    footnote(general = c("The model estimation included a local linear trend, triangular Kernel weights, and a data-driven bandwidth selection (mean squared error optimal bandwidth). The table displays the local average treatment effect in percentage points for individuals with a standardized blood pressure of zero and CIs resulting from robust bias-corrected standard errors clustered at the individual-level. The bandwidth was varied by +/- 0.2 in intervals of 0.05, column (5) contains the results."," Abbreviation: CI = confidence interval"),general_title="",threeparttable = T, fixed_small_size = T)%>%
    save_kable(filename, format = "latex")

```